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1.  RAY  TRACING  IN  THE  ARCTIC  OCEAN 


Ray  Tracing  Program 

In  the  Arctic  Ocean  it  is  known  that  below  ,the  upper  layer  of 
several  hundred  feet  the  velocity  gradient  is  fairly  uniform  and 
moderate  with  respect  to  depth  and  with  the  exception  of  extreme 
low  frequency  propagation  ray  tracing  affords  a  very  convenient 
and  effective  technique  for  calculating  the  acoustic  intensities. 

Aside  from  a  few  minor  details,  a  ray  tracing  program  is  now 
completed  and  is  capable  of  performing  the  following  tasks. 

Given  the  locations  of  a  number  of  sound  sources  (longitude, 
latitude,  and  depth)  and  the  locations  of  a  number  of  receiver 
stations  (longitude,  latitude,  and  depth),  the  program  evaluates 
the  acoustic  intensities  at  each  receiver  station  corresponding 
to  each  of  the  given  sources  and  the  topography  of  the  Arctic 
Ocean.  Input  data  consists  of  a  typical  temperature  and  salinity 
profile  and  the  RMS  ice  roughness  over  the  Arctic  Ocean  for 
computing  the  surface  reflection  losses. 

A  typical  ray  diagram  and  velocity  profile  using  recent  temper¬ 
ature  and  salinity  data  at  Base  Camp  (72-0. 20N,  148-35. 2W)  are 
plotted  and  included  in  this  report.  The  units  of  depth  ar.d 
range  are  in  feet  and  velocity  in  feet  per  second.  The  source 
is  1000  feet  below  the  surface  and  the  grazing  angles  of  the  rays 
have  one  degree  increment. 


Velocity  Profile  from  Temperature  and  Salinity  Data 

While  the  most  widely  used  formula  for  computing  the  velocity 
from  temperature  and  salinity  is  the  Wilson's  ("Equations  for 
the  Speed  of  Sound  in  Sea  Water",  JASA,  Vol.  32,  No.  10,  October 
1960,  page  1357),  there  is  a  shorter  and  simpler  formula  due  to 


Leroy  ("Development  of  Simple  Equations  for  Accurate  and  More 
Realistic  Calculations  of  the  Speed  of  Sound  in  Sea  Water," 

NATO  TR  No.  128  Salant  ASW  Research  Centre,  LaSpezia,  Italy, 
November  1968).  In  this  shorter  formula  there  is  no  need  to 
carry  out  a  separate  calculation  for  underwater  pressure  as  is 
necessary  when  using  Wilson's  (Leroy:  "Formulas  for  the  Calcu¬ 
lation  of  Underwater  Pressure  in  Acoustics,"  JASA  Vol.  44,  No.  2 
August  1968,  page  65). 


Both  Wilson's  and  Leroy's  formulas  have  been  used  to  make  some 
velocity  comparison  studies  from  typical  temperature  and  salinity 
profiles.  The  results  seem  to  indicate  that  the  difference  in 
velocity  is  quite  small,  being  of  the  order  of  one-tenth  of  one 
per  cent  at  large  depth.  In  the  ray  tracing  program,  a  choice 
can  be  made  of  either  Wilson's  or  Leroy's. 


Arctic  Topography 

Approximately  5000  data  cells  in  the  program  are  used  to  describe 
the  depth  of  the  Arctic  Ocean  bottom.  This  covers  the  rectangular 
area:  the  midpoints  of  the  four  sides  have  the  following  lati¬ 

tudes  and  longitudes:  70°N  0°E,  75°N  90°E,  70°N  180°E,  and  75°N 
90°WU  Each  data  cell  represents  the  average  depth  of  the  ocean 
bottom  for  an  area  of  approximately  35  miles  square. 

The  data  are  based  on  the  Arctic  Bathymetry  map  as  compiled  by 
the  Canadian  Hydrographic  Service,  1966-1967  (72°N  0°-90°W  and 
72  N  90°-180°W)  and  the  North  Polar  Chart  compiled  by  the  Hydro- 
graphic  Department  of  the  British  Navy,  1969. 

A  barthometic  computer  plot  showing  the  approximate  depth  of  the 
Arctic  Ocean  bottom  is  included  in  this  report. 
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2.  RELATED  THEORETICAL  DEVELOPMENT 
Corrections  to  Ray  Acoustics 

Let  (xi,  y i,  Zj)  be  rectangular  Cartesian  coordinates  and 
t  be  the  time.  Let  £  be  some  characteristic  length  (e.g., 
water  depth),  and  define  the  non-dimensional  coord,  nates 
by  x  =  xi /£,  etc.  Consider  the  wave  equation 

£2 

s  p-  *tt  (2.1) 

where  <p  is  the  velocity  potential,  c  =  c(x,y,z),  a  function 
of  the  space  coordinates.  Let  c0  be  some  constant  which  can 
be  taken  approximately  as  the  average  value  of  c  over  the 
region  of  interest. 


For  periodic  motion  we  assume 


<t>  -  Re{u  elwt }  (2.2) 

where  w  is  the  (constant)  angular  frequency;  u(x,y,z)  is  the 
complex  amplitude.  Substitution  of  (2.2)  into  (2.1)  shows 
that  u  satisfies  the  equation 


where 


Au  +  k2n2u  *  0 


(2.3) 


A  « 


a2  a2  92 

8X2  +  9y2  +  3Z2 


k  =  is. 

co 


and 


n 


1° 

c 


The  quantity  k  is  non-dimensional,  and  in  practice  k2  is 
usually  >>L  For  example,  at  a  frequency  of  20  Hz  and  with 
a  depth  of  2000  feet,  k2  is  approximately  2500. 


In  acoustic  propagation  problems  one  is  often  dealing  with  a 
sequence  of  progressing  waves  of  approximately  constant  phase; 
for  example,  plane  waves,  u  *■  e^x,and  spherical  waves,  u  se’^r/r* 
For  such  problems  a  natural  substitution  is 

u  -  a  (x ,  y,  z)eik0(x'y'z)  (2.4) 
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where  the  function  8(x,y,z)  is  to  provide  the  travelling 
wave  character.  Substitution  of  (2.4)  into  (2.3)  gives 

k2 [n2- (gx2+3y2+gz2)] a  +  ik [2 (axBx+ay8y+az&z)  +  aAg]  +  Aa  =  0 

(2.5) 

For  large  k  the  first  term  is  dominant. so  that  as  a  first 
approximation 


ex2  +  ey2  +  3z2  =  n2  0.6) 

This  is  the  eikonal  equation,  and  its  characteristics  are 
the  acoustic  rays.  If  s  is  the  distance  along  a  charac¬ 
teristic,  the  differential  equations  of  a  characteristic 
be  ~-ome 


dx  _  3x_  dy_  _ 

ds  n  ’  ds  n  *  ds  n 


d3x 

dS  =  "X* 


d3y 

ds  "  ny' 


dBz 

y- —  s  n 
ds  z 


dB  _ 
dF  58 


n 


(2.7) 


from  which  the  conventional  ray  equation  follows: 


d  2x 
ds  2 


1  dx  f  dx 
n  ds  ^Xds 


+  n. 


d£ 

ds 


+  n, 


dz> 
ds  * 


1 

+  — 
n 


n. 


(2.8) 

The  term  that  is  next  important  in  (2.5)  can  also  be  put 
to  zero  to  get  a  second  approximation  by  choosing  a  such 
that 

2(ox8x  +  ayBy  +  az8z)  +  aAB  ■  0  (2.9) 
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mm* 


From  (2.7)  and  (2.9)  it  foil 


ows  along  a  ray 


2n^  +  ctA3  =  0 
ds 

i,e,»  div(a2  grad  3)  =  0 

and  by  the  divergence  theorem  for  a  closed  surface  S 

/  a2  grad  3  •  dS  =  0  (2.12) 

s 

Consider  a  section  of  a  ray  tube  and  make  use  of  (2.12).  The 
contribution  to  the  integral  from  the  tube  sides  is  zero  since 
grad  3  is  in  the  direction  of  the  ray;  so  that  if  the  two  end 
sections  of  the  ray  tube  be  given  the  indices  1  and  2,  we  have 


(2.10) 

(2.11) 


/  a2ndS=  J  a2ndS 
S1  S2 


(2.13) 


Since  the  rate  of  the  energy  propagation  per  unit  area  is  propor 
tional  to  (amplitude) 2/c ,  (2.13)  implies  that  the  rate  of  energy 
propagation  is  constant  along  a  ray  tube. 


This  leads  to  the  conclusion  that  the  use  of  ray  theory  is  equiv¬ 
alent  to  the  omission  of  the  last  term  in  (2.5).  One  way  of 
examining  the  adequacy  of  the  ray  theory  is  to  solve  the  problem 
by  ray  theory  to  obtain  o  and  3  at  each  point  and  to  examine  the 
magnitude  of  Aa.  Another  approach  is  to  iterate,  compute  Act  =  f, 
say,  and  recompute  a  by 

2nj|  +  aA3  +  f  =  0  (2.14) 


A  comparison  of  the  revised  values  of  a  should  provide  a  criterion  for 
the  adequacy  of  the  ray  approximation. 


Such  iterations  can  be  carried  out  more  systematically  by  usin^  the 
asymptotic  expansion  below. 
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«  =  a(0)  +  i_  aO)  + 

i  k 


&  “(2)* 


(2.15) 


Substitution  into  (2.5)  with  5  still  satisfying  (2.6)  yields  the 
sequence  of  equations: 


2n 


da 


(0) 


ds 


+  a^Ae  = 


2n 


da 


(1) 


dF 


+  a^Afi 


+  Aa 


(0) 


2n 


da 


(2) 


JT 


♦  a  AS 


+  Aa 


(1) 


(2.17) 


so  that  ....  may  be  determined  in  succession.  Note 

that  the  ray  paths  are  unaffected  by  this  iterative  procedure. 
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3.  SUBCONTRACT  WORK  TO  LAMONT-DOHERTY 

The  principal  effort  at  Lamont-Doherty  has  been  to  develop  a 
rapid,  accurate  method  of  computing  propagation  loss  as  a  function 
of  range  in  the  ice  covered  Arctic  Ocean.  Important  input  param¬ 
eters  to  the  propagation  model  are  ice  roughness,  bottom  topography, 
and  the  velocity  structure  as  a  function  of  depth  in  the  ice,  ater, 
and  bottom  sediments.  Computational  speed  is  of  the  utmost  impor¬ 
tance  to  evaluate  variations  of  these  parameters  on  propagation 
loss  as  a  function  of  range  since  a  number  of  models  must  be  con¬ 
sidered  . 

Such  a  rapid  computational  method  is  the  Fast  Field  Program  (FFP) 
technique  discovered  by  Marsh  (1967).  The  FFP  was  implemented 
for  computation  on  a  digital  computer  by  Di  Napoli  (1971)  for 
propagation  in  an  all  liquid  wave  guide.  In  the  Arctic  Ocean  solid 
layers  as  well  as  liquid  layers  must  be  considered,  and  a  conven¬ 
ient  starting  point  for  the  FFP  is  the  integral  solution  derived 
by  matrix  methods  by  Kutschale  (1970,  1972)  for  propagation  from 
point  harmonic  sources  in  a  multilayered  liquid-solid  half  space. 

In  the  FFP  technique  the  integral  solution  is  evaluated  rapidly 
as  a  function  of  range  by  numerical  integration  employing  the 
Fast  Fourier  Transform  (FFT) .  Singularities  in  the  integrand 
corresponding  to  the  normal  mode  poles  are  remove-  from  the  axis 
of  integration  by  including  an  attenuation  factor  in  each  layer. 

The  attenuation  in  each  layer  can  represent  the  effects  on  propa¬ 
gation  loss  as  a  function  of  range  of  absorption  of  sound  in  the 
water,  attenuation  by  the  rough  ice  boundaries,  and  attenuation 
of  sound  by  the  bottom. 

Two  computer  programs  were  written  in  Fortran  IV  to  evaluate  the 
liquid-  and  solid-bottom  integral  solutions  derived  by  Kutschale 
(1970,  1972).  A  comparison  of  the  FFP  computations  with  those 
computations  by  normal  mode  theory  from  the  corresponding  integral 
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solutions  are  nearly  identical,  but  the  FFP  technique  is  far  more 
convenient  and  at  least  an  order  of  magnitude  faster  since  the 
computations  are  done  directly  from  the  integral  solution  without 
first  computing  the  roots  of  the  frequency  (period)  equation  and 
then  summing  the  normal  modes. 


Computations  done  to  date  by  the  FFP  are  in  close  agreement  with 
measurements  of  propagation  loss  in  the  central  Arctic  Ocean. 
Present  effort  is  a  further  comparison  of  propagation  loss  compu¬ 
tation.1;  against  experimental  data  and  a  comparison  between  FFP 
and  computations  by  ray,  normal  mode,  and  the  WKB  mode  theory. 

An  evaluation  is  in  progress  of  the  effects  on  propagation  loss 
as  a  function  of  range  of  thermal  and  salinity  gradients,  surface 
and  bottom  reflections,  and  bottom  topography. 
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